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Abstract 

The paper provides full algorithmic details on switching to the continuation of 
all possible codim 1 cycle bifurcations from generic codim 2 equilibrium bifurcation 
points in n-dimensional ODEs. We discuss the implementation and the performance 
of the algorithm in several examples, including an extended Lorenz-84 model and a 
laser system. 



1 Introduction 



Consider a system of differential equations depending on two parameters 

x = f(x, a), (x, a) G R n x R 2 , (1) 

where / is smooth. In general, there are bifurcation curves in the a-plane, 
at which the system exhibits codim 1 bifurcations, for example, fold or Hopf 
bifurcations of equilibrium points. Moreover, genericaliy, one expects points 
of codim 2 bifurcations, where several curves corresponding to codim 1 bifur- 
cations intersect transversally or tangentially. A codim 2 point is of particular 
interest if it is not only the origin of some equilibrium bifurcation curves but 
also of some curves corresponding to bifurcations of periodic orbits (cycles). 
Such points can be detected by purely local analysis of equilibria and then 
be used to establish the existence of limit cycle bifurcations and other global 
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phenomena that could hardly be proved otherwise. That is why codim 2 points 
are often called the "organizing centers" in applied literature. 

The theory of codim 2 bifurcations of equilibria in generic systems (CQ) is well- 
developed (see, for example, [I], [12], |18j). There are five well-known codim 
2 equilibrium bifurcations: cusp (CP), Bautin (generalized Hopf, GH), double 
zero (Bodanov-Takens, BT), zero-Hopf (ZH), and double Hopf (HH). It follows 
from their analysis that branches of nonhyperbolic limit cycles can emanate 
from GH, ZH, and HH points only. More precisely, a codim 1 bifurcation curve 
LPC, along which a cycle with a nontrivial multiplier /ii = 1 exists, emanates 
from a generic GH point, while codim 1 bifurcation curves NS, along which 
a cycle with a pair of multipliers = e ±%e exists, are rooted at generic ZH 
and HH points. Notice that NS is used to denote both Neimark-S acker and 
neutral saddle cycles where ^1^2 = 1 and that no period-doubling curves can 
emanate from generic codim 2 equilibrium bifurcations. 

Obviously, the application of these theoretical results to realistic models ([I]) 
is impossible without numerical tools. The numerical analysis of a codim 2 
equilibrium bifurcation includes: 

• detection and location of the point in a branch of a codim 1 bifurcation; 

• computation of the coefficients of the normal form of the restriction of ([T|) to 
the critical center manifold at the bifurcation parameter values and checking 
the nondegeneracy conditions; 

• verification of the transversality of the given family (CQ) to the codim 2 bi- 
furcation manifold and establishing a correspondence between the unfolding 
parameters of the normal form and original system parameters a; 

• computing accurate approximations of the codim 1 curves in the a-space 
and the corresponding singular orbits in the x-space near the bifurcation, 
sufficient to initialize the numerical continuation of these codim 1 curves 
using only local information available at the codim 2 point. 

While the first two problems were studied in detail (see, [3] and references 
therein) and have been implemented into the standard bifurcation software 
content [TU] and matcont [5J, two last issues received much less attention 
in the numerical analysis literature, even if bifurcations of nonhyperbolic cycles 
are concerned. The present paper is aimed at bridging this gap by providing 
full algorithmic details on switching to all possible codim 1 cycle bifurcations 
from generic GH, ZH, and HH codim 2 points. 

One way to set up a computational switching procedure is to consider a smooth 
normal form for the codim 2 bifurcation including the parameters (3 G M. 2 

w = G(w, /?), G:f" c xl 2 ^ R nc . (2) 

For all codim 2 equilibrium bifurcations these normal forms are known. Sup- 
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pose that an exact or approximate formula is available that gives the emanat- 
ing codim 1 bifurcations for the normal form (j2j). In order to transfer this to 
the original equation (JTJ) we need a relation 

a = V(J3), V : R 2 -> R 2 (3) 

between the unfolding parameters (3 and the given parameters a. In our con- 
text, V will be linearly approximated. Moreover, we need a center manifold 
parametrization 

x = H(w, (3), H : R" c xl 2 ^ R", (4) 

that incorporates (3. Taking (j3J) and (jlj) together as (x,a) = (H(w, 0), V(f3)) 
yields a center manifold for the suspended system x = f(x,a),a = 0. The 
invariance condition for the center manifold now turns into a homological 
equation: 

H w (w, (3)G(w, (3) = f(H(w, (3), V((3)), (5) 

which we can solve by a recursive procedure based on Fredholm's solvability 
condition that will give the Taylor coefficients of G and H with respect to w 
and (3. We assume the Taylor series of G to be known as 

G(w,(3)= £ -L^^m 
M+H>i ZA/i ' 

and the Taylor series of H and V to be unknown 

H(w,(3)= -^-jKrvrp, V({3)=J2-, V ^- 

Here v and \x are multi-indices. For \x = this reproduces the critical normal 
form coefficients first computed in [T7j, while the coefficients with \fi\ > 1 yield 
the necessary data on the parameter dependence. 

To summarize, a bifurcation point is detected within a certain small tolerance. 
As the prediction depends on the initial point, this translates into small errors 
of the predicted curve. If we start close enough to the actual new curve, any 
point will converge to it and in general one expects a convergence cone |15j . 
If we parametrize the predicted curve by e, the initial amplitude e is to be 
chosen to be within the convergence cone, see also Figured! 

This procedure is adopted from [3], where it has been applied to the derivation 
of the asymptotics of the fold and Hopf curves rooted at CP and BT codim 
2 points, as well as that for a homoclinic orbit to a saddle emanating at the 
BT-point. Recently, this technique has been successfully used for switching at 
codim 2 fixed points of maps to the continuation of nonhyperbolic periodic 
orbits rooted there [8]. Note that a similar procedure was suggested in [13] . 
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Fig. 1. Sketch of the switch in the case of a GH bifurcation. A predicted point along 
P (e) (a) in the cone will converge to the LPC-curve, outside (b) it will not. 

without using the Fredholm condition, and carried through in the ZH-case 
in [14] . where, however, no asymptotics of codim 1 curves were derived. Fi- 
nally, we point out that the problem of switching to the LPC-curve at the GH 
bifurcation has been briefly discussed in [3] in a setting without the Taylor 
expansion in (3. 

The paper is organized as follows. In Section [2] we revise smooth parameter- 
dependent normal forms on center manifolds for the considered codim 2 bifur- 
cations, i.e. give G(w,P) in GH, ZH, and HH cases, and give the asymptotic 
expressions of the branches of nonhyperbolic cycles in these normal forms. 
Then we perform the described above reduction procedure and derive the nec- 
essary coefficients g u ^,h Ufl , and t> M in terms of F and its derivatives. These 
coefficients are finally used to set up predictors for these branches in the orig- 
inal system (CQ). An implementation of the resulting formulas in the software 
matcont is discussed at the end of this section. Section [3] presents several 
applications of the developed technique to known ODE models, an extension 
of the Lorenz-84 system and a laser model, where we compare the asymptotic 
formulas for the cycle bifurcations with numerically computed LPC- and NS- 
branches. A discussion of existing results and open problems in switching to 
homoclinic branches at ZH and HH bifurcations is given in Section HI 



2 Asymptotics and the Center Manifold 

2.1 The 'new' curves 

The parameter-dependent normal forms are known and can be found in the 
standard texts, e.g. [18J. As the normal form and the asymptotic expressions 
are the necessary theoretical ingredient, we present these here. 
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2.1.1 Generalized Hopf 



Near a GH bifurcation the vector field restricted to the center manifold is given 
by 

w = \(P)w + Cl (P)w\w\ 2 + c 2 ((3 2 )w\w\ 4 + 0(\w\ 6 ), u>GC, (6) 

where A(0) = iu, and this bifurcation is characterized by d\ = 3?(ci(0)) = 
and d 2 = dt(c 2 {0)) 7^ 0. A curve LPC of fold bifurcation of limit cycles emanates 
from this point. Let us write w = pe 1 ^, A(/3) = iu + (3\ + ib\((3) + 0(\fl\ 2 ) and 
dt(ci(P)) = 02 + 0(\/3\ 2 ). If we now truncate the normal form to fifth order in 
w, then the curve LPC is given by 



£,A = d 2 e\f3 2 = -2d 2 e 2 . (7) 



2.1.2 Zero-Hopf 



Near a ZH bifurcation the vector field restricted to the center manifold is given 
by 









-( 







A + /200a; 2 + /011M 2 + /300a; 3 + /m^kl 2 

{iu(/3) + f3 2 )w + ^uoanw + g 2 iox 2 w + g 02 iw\w\ 2 



where (x, w) G M x C. An extra Neimark-Sacker (torus) bifurcation of limit 
cycles (NS) occurs if 3?(<7no)/oii < 0. 



The asymptotic expression is 



/ill + 2(7021 _2 n f 2 



J 111 1 -JU^i An £ V 

p = £) x = — £ , ft = -/one 

^/200 

„ 2(g(gno) - /200Mff02l) + %(fliio)/m 2 

P2 ^7 e • 

^J200 



(9) 



This agrees with a formula given in [TJ. 



2 J. 3 Double-Hopf 



For a HH bifurcation the dynamics on the center manifold is governed by the 
following normal form: 



' Wi \ I {uJi{0) + Pi)wi + f2iooWi\wi\ 2 + /1011W1IW2I 2 



-o(\\( Wl ,w 2 w 



{w 2 J \(iu 2 (P) + h)w 2 + gui W2\w 1 \ 2 + g 00 2iw 2 \w 2 \ 

(10) 

where (w±, w 2 ) G C x C. Then there are generically two half- lines along which 
there is a NS bifurcation of limit cycles. In polar coordinates w\ = pie 1 ^ 1 , w 2 = 
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p2^ 2 their asymptotics are given as 



(pi, p 2 , A, £2) = (e, 0, -3?(/ 2 ioo)£ 2 , -^(^1110)^) , 
(pi, pa, A, /%) = (0, e, -3?(/ion)^ 2 , -&(0oo2i)e 2 ) • 



(11) 
(12) 



2.2 Coefficients of parameter- dependent center manifolds 



We assume that the critical normal form coefficients are known (see [T7] and 
[3] ) and give here only parameter- related coefficients h vll from the homological 
equation. These provide in each linear approximation to the parameter 

transformation . 



2.2.1 Generalized Hop f 

Here we closely follow the idea outlined in [3] . We first expand the eigenvalue 
and the first Lyapunov coefficient in the original parameters a and collect 
the equations to obtain the transformation to the unfolding parameters (3. 
Alternatively, one can normalize already in ({TBI to obtain an orthogonal frame 
from these equations and obtain scalings from the higher order equations (|14p . 
Below we have p = (10), (01) as indices and viq = (1, 0), t>oi = (0, 1) as vectors. 

The first two equations (actually four) coming from ((Sj) are 



The first equation is nonsingular and from the second we find 71.^ using the 



Ah, 



OOfj, — 



(13) 



(A - iujI n )h WlJi = j ltfi q - A t (q, u M ) - B(q, h 00)M ) 
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Fredholm alternative. The other systems from (j3j) are 



(A - 2iuI n )h 20/Ji 
(A - iul n )h 21fl 



2/i2ooo7i,^ - [C(q, q, h 00fJi ) + 2B(q, h 10li ) + B(h 200 o, h 00ll ) 
+-Bi(g,g,u M ) + Ai(/i 20 oo,^)] , 
23?(7i )At )/i U oo - [C(g, g, /iqo^) + £(Anoo, h 00fl ) 
+B(q, h Wli ) + ^oi/i) + 5 i(<?, 9> *V) + ^1(^1100, *v)l , 
272,^g + ^2100(271,^ + 7i,m) + 2^10/^1 
- [D(q, q, q, h 00fl ) + 2C(q, h noo , h 00ll ) + 2C(q, q, h 10ft ) 
+C(q, q, hoin) + C(h 20Q0 , q, h m ^) + 2B(q, h llfl ) 
+2B{hixoQ, hion) + B(h 2 ooo, h mtJ ) + ^(^2100, ^00/J 
+B(h 20 ^,q) +Ci(q,q, q, v M ) + 2B X (h noo , g, 

+B-L (/l 2 000, q, Vfi) + Ai (/l2100, U/i)] , 

(14) 

The first two are nonsingular and with the Fredholm alternative we find 72.^. 
The parameter transformation ([3]) is given by 



a 



ft 



^7i,io Ti.oi 

^72,10 72,01 



-1 



(3. 



(15) 



2.2.2 Zero-Hopf 

This case is also treated in [J3], however with only one parameter and for 
hyperbolic periodic orbits. Thus our computational scheme is different. We 
list only the necessary equations. 

(a) A[h 00010 , /iooooi] = [qi, 0]-Ji[vio, v Q1 ], 

(b) A[h W010 , h W001 ] = [h 20000 , 0] - Ai(q u [v 10 , v 01 )) 

-B(q u [/ioooio, h 00001 ]) (16) 

(c) (A - iul n )[h 0ww , V001] = [^11000, 92] - 4i(?2, [vio, voi}) 

—B(q 2 , [/ioooiO; ^00001]) 

In contrast to the other cases, here the first system is already singular. Taking 
the inner-product with the adjoint null-vector we obtain the new orthogonal 
frame 



7 = (7i,72) =pJJi, 



7/II7II > s 2 = (-72,7l), 



Via = si + 5is 2 ,^oi = S 2 s 2 . 



(17) 
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Polynomial terms in the normal form (jSj) like j3\x are also resonant, but they 
can be eliminated by hypernormalization. After solving (fISI a) with a bordered 
matrix, see [IT] , still a multiple of q\ may be added to /^oooio- We use this to 
perform hypernormalization. Let us write 



n = -A 



i V1 I <?i ~ J i s i 




/ 



r 2 = -A 



INV 



-JlS2 





where A INV indicates the use of the bordered matrix, then we can write 



'00010 



00001 



S 2 r 2 + 5 A qi, 



for some <5's. Then by applying the Fredholm alternative to (1161 b. c) we can 
solve for all 5's at once. 



LL 



U(LL) 



Si 



\$3 



V 





(pi, Ai{qi,n) + B(q 1 ,r 1 )) 
(p 2 ,A 1 (q 2 ,r 1 ) + B(q 2 ,r 1 )) 



where 



LL 



(p 1 ,A 1 (q 1 ,r 2 ) + B(q 1 ,r 2 )) 2f< 



200 



v (P2,A 1 (q 2 ,r 2 ) + B(q 2 ,r 2 )) g m 



2.2.3 Double Hopf 

Although high- dimensional, this case can be treated in a relatively simple 
manner. Using the same notation as for the generalized Hopf from §5§ we get 

^OOOOa* = —Jl v /ij 

(A - zo;i4)/iiooo^ = 7i,m<?i - A(<?i,^ M ) - B(q l7 h 0000 ^), (19) 
(A - iuJ 2 I n )h 0010f , = 7 2>Ai g 2 - A(<?2, v^) - B(q 2 , /i ooo M )- 

As the first equation is non-singular, formal substitution of /iooooio and /loooooi 
and the Fredholm alternative leads to the same transformation (fT5|) from 
unfolding to the system parameters. 
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2.3 Implementation of the Predictors 

We have implemented our switching routines in matcont [6] . For the contin- 
uation of LPC and NS curves it uses a minimally augmented defining system 
|lbj . i.e. we need to supply an approximation of the limit cycle, the period and 
the parameters. The parameters follow from applying the inverse transforma- 
tion to ([3]). There is always one dynamic variable ip giving a free phase shift 
along the bifurcating limit cycle with a period ^hj- For the initial cycle we 
make an equidistant mesh ip = 2nn/N, n = . . . N where N + 1 is the number 
of mesh points. Let q denote the eigenvector corresponding to the eigenvalue 
iui, then points on the limit cycle are given by xo + e(qe 1 ^ + qe~ 1 ^). Similarly, 
terms as e 2 h 2 oe 21 ^ and £ 2 /iooio are included. An internal routine of matcont 
then adapts this limit cycle on an equidistant mesh to a mesh defined at the 
non-equidistant collocation points. 

For the NS curves the system is augmented with the real part k of the multi- 
plier. In this case the normal forms (IB]) . (fTOl) also define a second rotation with 



matcont uses Moore-Penrose continuation for which also a tangent vector 
to the bifurcation curve is needed. This tangent vector is easily obtained by 
differentiating the predictor w.r.t. e. 

Below we list some case-specific details. 

2.3.1 Generalized Hop f 

The period is given by T = 2ir/u + {2d 2 b h2 - 9f(ci (O)))^ 2 , with b lj2 = f|. The 
parameters are given by a = a + V(0, — 2d 2 e 2 ) T . 

Note that for a ^-approximation also seventh order derivatives would be 
needed; this follows from Remark 3.3.2 in [2T] . Therefor we restrict to 0(e 3 ) 
in the implementation. 

2.3.2 Zero-Hopf 

In the continuation we also need to provide the period and the multiplier. 
Approximating formulas are defined as follows where x,(5i,f3 2 are as in ([9]) 




T = 2n/uo(0) - (wift + u 2 f3 2 + $f(jg no )x) - %(g 2i)e 2 , 
k = l-(4n^(g no )f 011 )(e/u; ) 2 . 



(20) 
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2.3.3 Double Hopf 

Approximating formulas for the period and the multiplier on one branch are 
given by 

T = ^kT^ fc = cos(T(^ 2 + ^ 2 )), 

(dux, dcu 2 ) = -S(7i7 2 ) T (K(7i7 2 ) T )~ 1 K(/2ioo, 9nio) T + 3(/ 2 ioo, 9mo)- 
and similarly for the other branch. 



3 Examples 

3.1 New curves in an extension of the Lorenz-84 model 

The first example is an extended version of the Lorenz-84 model. A bifurcation 
analysis of this model was presented in [22123] • In this system X models the 
intensity of a baroclinic wave and Y and Z the sine and cosine coefficients of 
the wave. This model may be extended with a variable U to study the influence 
of external parameters such as temperature and the model then shows several 
limit cycle bifurcations [20J. It has the form: 

' X = —Y 2 - Z 2 -aX + aF- 7 f/ 2 

Y = XY - I3XZ -Y + G 
< . (22) 
Z = (3XY + XZ — Z 

^ JJ = -5U + >yUX + T 

The parameters F and T are varied while we fix a = .25, (3 = 1, G = .25, 5 = 
1.04,7 = .987. The bifurcation diagram displays one fold bifurcation and two 
Hopf bifurcation curves, see Figure [21 We find all codim 2 points of equilibria, 
in particular GH, ZH and HH. 

We have applied our switching routines to all three emanating curves, since 
the iVS* bifurcation from ZH is a neutral saddle. The predictions in parame- 
terspace are shown in Figure [2] next to the numerically continued curves. The 
predicted points were used as a starting point point for the continuation of 
these limit cycle bifurcations, which shows that our approach works. Another 
numerical check is provided by inspecting the tangent vector, which we pro- 
vide together with a first point. When we find a second point on the curve 
by continuation and adapt the defining system, we will obtain a more precise 
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Label 


XT' 

r 


1 


Normal Form coefficients 


GH 


2 3763601 


050197432 


do — 1558012 


HH 


2.5332211 


.026273943 


P11P22 = -1,0 = -3.648550, 5 = -1.052987 








G = 1230.630, A = -210.861 


ZH 


1.2834193 


.000126541 


8 = 1,9 = 0.3715145, E = -1 



Table 1 



Parameter values of F and T at the bifurcation points in Figure [2] together with 
normal coefficients (scaled, see [18]). 




Fig. 2. Bifurcation diagram of the Extended Lorenz-84 model. Symbols denote LP 
Limit Point, H Hopf, LPC Limit Point of Cycles, NS Neimark-S acker, GH Gen- 
eralized Hopf, HH Double Hopf, ZH Zero-Hopf, £>T=Bogdanov-Takens. Dashed 
lines show the predicted new curves; (a) Zoom near the GH point, (b) Zoom near 
the HH point. 

tangent vector. For a small continuation step, this tangent vector and the pre- 
dicted one should be close. Indeed, for the examples reported here, the first 
digits always coincided. 

Finally we present some measure of the error of the switching routines as a 
function of the initial amplitude e, see Figure [3] and its caption. Interestingly 
this Figure represents the idea of Figure [TJ Using a small initial amplitude e 
may not work due to a numerical error in the calculated codim 2 point, on 
the other hand e must not be taken too large for the approximation to remain 
valid. 
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Fig. 3. Error measures: (a) The residual R of the first Newton-step, (b) The distance 
between the predicted and the first corrected point. For this we left out the tangent 
vector in the continuation, matcont then tries to correct the first point immediately 
instead of starting the continuation. We have taken 20 mesh and 4 collocation points 
and e € [10 , .2]. Data shown if predictor converged. 

3.2 Switching in a Laser model 



In [21] a single-mode inversionless laser with a three-level phaser was studied 
and shown to operate in various modes. These modes are "off" (non-lasing) , 
continuous waves, periodic, quasi-periodic and chaotic lasing. In particular 
the boundary of the region of chaos seems to be defined by several limit cycle 
bifurcations born from several codim 2 equilibrium bifurcations. Thus we want 
to start with our routines such boundary computations without first doing 
simulations and limit cycle continuations in this 9-dimensional system. 

The model is 9-dimensional system given by 3 real and 3 complex equations: 

Paa = Ra~ ffllfaab - <b) + ^p(^ac ~ <J), 

f) bb = R b + ffil((T ab - <7* b ), 
< [26) 

O-ab = -(7l + «A/)(7 a6 - ^{Vtl{p aa ~ Pbb) ~ ^pO"cfe), 
0~ac = -(72 + i\)0- ac ~ I {?L p {2paa + Pbb ~ 1) ~ ^cb) , 
t CTcfe = -(73 + ~ A p ))(T c6 - \ (£Va C - QpO-ab), 

with R a = -.505p aa - A05p bb + .45, R b = .0495p aa - .0505p 66 + .0055 and A, := 
A cav + g^(o- ab )/Qi. The parameters are fixed at 71 = .05, 72 = .25525, 73 = 
.25025, 7 ccro = .03, g = 100, A p = while A cav and Q p are varied to study 
several detuning effects. For more details, see [24"] . 

We have reproduced a part from the bifurcation diagram which corresponds 
to continuous wave and periodically pulsating solutions, i.e. with Qi 7^ 0, see 
Figure HI As the system has Z 2 -symmetry the same bifurcations are found for 
A cav — * —A cav . For clarity of the figure we do not display these here. We list 
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Fig. 4. Bifurcation diagram of the inversionless laser. Hopf curves (denoted by H) 
are dotted, Limit cycle bifurcations are denoted by LPC Limit Point of Cycles, NS 
Neimark-S acker. Dashed lines show the predicted new curves. 

the codim 2 points in Table I3T21 The normal form coefficients of HHi confirm 
the claim of [21] that the most complicated type was encountered; only the 
3-torus is (un) stable. This is also confirmed when we continue the Neimark- 
Sacker bifurcations. For HH 2 the NS curves are not in the same quadrant 
defined by the Hopf curves, while they are for HHi. All cycle bifurcations 
where computed with 20 mesh points and 4 collocation points and the initial 
amplitude was set to e — .001, which worked immediately in all cases. Let us 
remark that one LPC curve connects GH2 and GHz points and stays close to 
the Hopf curve. Similarly, a NS curve starts at HHi, becomes neutral between 
two 1:2 resonances and ends at HH 2 . It would have taken much more effort to 
find this feature otherwise. 



4 Discussion 



This paper contributes to the bifurcation analysis of codim 2 singularities 
of equilibria in multidimensional ODEs by providing explicit predictors for 
branches of nonhyperbolic cycles emanating from these bifurcations. We have 
tested it on several examples with good results. We believe that this work will 
further facilitate automated analysis of nonlinear systems. However, we like 
to mention that we also tried the double Hopf point in a model for the lateral 
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-Laoei 




A 

'-^cav 


Normal Form coefficients 


Grill 


*7 OOQ01 n 

i .zzooiy 


0.011400 


0,2 = — 4b.4yo0z 


LrU2 


r HOI K*7A 

O.UzlO I 4 


1.44DOO I 


J Q QI QI on 

U2 = o.ololoz 


LrH 3 


A QOAPidP. 

4.oz4U00 


l.Uoyoo i 


j i net 1 1 1 n 
$2 = lyo.lliy 




o qi oi on 


— o.zY oobo 


02 = — b.4bo4bo 


HHi 


5 0872QQ 


— 1 2362053 


„,,„„„ - _i # — — 07194543 ^ — — 13Q1412 








9 = .9595389, A = -2602.275 




3.555848 


-1.983857 


P11P22 = 1,9 = -.1179924, 5 = -26.59452 








9 = -10.81042, A = -2713.608 



Table 2 

Parameter values of f2 p and A cav at the codim 2 points together with normal coef- 
ficients (scaled, see [18]). 




(a) 10 log( £ ) (b) 10 log( £ ) 



Fig. 5. Error measures along the eight curves: log of the distance between the 
predicted and the first corrected point versus 10 log(e). (a) Along the Neimark-S acker 
curves (b) Along the LPC curves. This Figure again resembles the idea of Figure [TJ 

pyloric neuron PfTTJ] . Although we were able to switch to one branch and 
continue it without any problem, the Jacobian of the denning system along 
the second branch was numerically singular. In this model with multiple time 
scales probably a special numerical scheme is necessary. 

It is well known that branches of orbits homoclinic to hyperbolic equilibria 
are also rooted at BT, ZH, and HH codim 2 bifurcation points. The BT case 
has been treated in [2] (see also [3], where the computational setting is most 
close to the present paper). The corresponding predictor for the homoclinic 
branch is implemented in matcont. The problem of providing predictors 
for homoclinic branches rooted at ZH and HH points is more challenging. 
Some important results in this direction are obtained in [4.7.5J , where the 
systems reduced to the center manifold were considered. However, a complete 
set of formulas suitable for switching to homoclinic curves in these cases is 
still not available. For instance, in the ZH case the normal form ([HD exhibits 
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homoclinic bifurcations of saddle-focus equilibria in the parameter plane along 
a bifurcation curve with the linear approximation 
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provided that 3?(fl , no)/on < an d ^(fl'iio)/2oo < 0. Application of the above 
reduction to the parameter-dependent center manifolds in the ZH case yields 
an approximation to the bifurcation curve in the parameter plane. Now the 
challenge is to construct a suitable initial solution in state space. On this work 
in progress will be reported elsewhere. 

Another direction for future research is a problem of switching to secondary 
cycle bifurcations at codim 2 bifurcations of cycles in ([T]). Here a generalization 
of the periodic normalization technique from |16j to critical codim 2 cases and 
its extension to parameter-dependent systems in the spirit of [S] are required. 
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